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Abstract. Over the past few years, we developed a mathematically rigorous 
method to study the dynamical processes associated to nonlinear Forchhcimer 
flows for slightly compressible fluids. We have proved the existence of a geomet- 
ric transformation which relates constant mean curvature surfaces and time- 
invariant pressure distribution graphs constrained by the Darcy-Forchheimer 
law. We therein established a direct relationship between the CMC graph equa- 
tion and a certain family of equations which we call 3-Forchheimer equations. 
The corresponding results, on fast flows and their geometric interpretation, can 
be used as analytical tools in evaluating important technological parameters 
in reservoir engineering. 



I. Background 

The present report uses methods of differential geometry and integrable systems 
in modehng nonhnear flows in porous media. The classical flow, as defined by 
Darcy's equation [13l [HI O [8] represents a linear relation between the gradient 
of pressure and velocity field. Darcy's formulation of the momentum equation of 
motion essentially simplifies some complicated hydrodynamics in porous media, by 
using a permeability tensor as the major characteristic of the media at a point 
(pore, block, layer, etc) in space. 

Most of the up-scaling techniques in heterogeneous reservoirs are based on the 
assumption that on each scaling level, the constitutive momentum relation between 
the gradient of pressure and velocity field remains invariant, while the only param- 
eter that is reevaluated is the permeability tensor. 

Philipp Forchheimer, one of the founders of groundwater hydrology at the be- 
ginning of the 20th century, observed that for high velocity rates, Darcy's law is 
no longer valid. On the other hand, on a fine level, it is well-known that even slow 
flows can deviate from the linear ones [14]. 

Some newer theories suggest that it is precisely the inertial term that causes the 
deviation from the linear Darcy law, but the nature of these inertial forces is not 
fully understood. As it was pointed by Bear, [5] (see also [24l [23l [17] ) : "Most of the 
experiments indicate that the actual turbulence occurs at a Reynolds number at 
least one order of magnitude higher than the Reynolds value at which the deviation 
from Darcy's law is observed" . In some of the experiments, a deviation from Darcy's 
law was observed for Reynolds number Re = qd/ fj, ~ 10 (where d is some length 
dimension of the porous matrix, q a specific discharge per unit cross-sectional area 
normal to the direction of the flow, and fi is viscosity of the liquid). In some recent 
works, such as [23], it was experimentally observed that Darcy's law is not verified 
even for Re ~ 1, for samples of the porous rocks containing fracture. Latest results 
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suggest that even some low velocity flows in a highly heterogeneous reservoir may 
deviate from Darcy's linearity principle. 

In reservoir engineering, there are three most popular non-linear approximations 
of the field data, establishing a formula for the pressure drop AP in terms of the 
production rate Q: 

- the "two term" law - AQ + BQ"^ = AP; 

- the "power" law - CQ" + aQ = AP, 1 n < 2; 

- the "three term" cubic law - AQ + BQ'^ + CQ^ = AP. 

All three equations were originally introduced by Forchheimer in works published 
at the beginning of the 20th century. 

Following [24j|4l[l][3] and references herein we have established three types of gen- 
eralized nonlinear Darcy law, with permeability tensor depending on the gradient of 
the pressure function. We will show that these constitutive equations correspond to 
a generalized Forchheimer equation. Under some assumption about fluids, the gen- 
eralized Darcy-Forchheimer equations enable a reduction of the system of equations 
that governs the flow - namely to one parabolic non-linear equation for the pressure 
function. This parabolic equation displays some similarity to the constant mean 
curvature (CMC-graph) equation for surfaces. In this work we use this similarity to 
find a constraint on the Forchheimer flows, such that the pressure function can be 
regarded as a surface with given constant mean curvature, after a certain geometric 
transformation. Conversely, a graph with prescribed constant mean curvature can 
be interpreted as a pressure distribution of the flow subjected to a nonlinear Darcy 
law, similar to the Forchheimer equation. This geometric interpretation provided 
us a with a simple algorithm to compute the productivity index of the well, in a 
structured heterogeneous reservoir. The productivity index (PI) of the well is one 
of the fundamental concepts in reservoir engineering, defined as the ratio between 
the production rate and the difference between the well pressure and the average 
pressure in the reservoir. The PI characterizes the well performance with respect 
to the geometry of the hydraulic system. It is shown that, for some specific con- 
dition that is reasonable to either impose or approximate, the productivity index 
of the well in the inhomogeneous reservoir can be computed using a solution of a 
corresponding CMC graph equation. 

II. Introduction. 

II. 1. General Forchheimer equations. Darcy's law (for viscous fluid laminar 
flows) assumes that the total discharge is equal to the product between the medium 
permeability, the cross-sectional area of the flow, and the pressure drop, divided by 
the dynamic viscosity. After dividing both sides of this equation by the area, one 
obtains another way to express Darcy's law: namely that the filtration velocity (or 
Darcy fiux) is proportional to the pressure gradient, through a certain permeability 
coefficient. Darcy's equation, the continuity equation and the equation of state 
serve as the framework to model processes in reservoirs [m [8]. For a slightly 
compressible fiuid, the original PDF system reduces to a scalar linear second order 
parabolic equation for the pressure only. The pressure function is a major feature 
of the oil or gas filtration in porous media, which is bounded by the well surface 
and the exterior reservoir boundary. Different boundary conditions on the well 
correspond to different regimes of production, while the condition on the exterior 
boundary models flux or absence of flux into the drainage area. All together, the 
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linear parabolic equation, boundary conditions and some assumptions or guesses 
about the initial pressure distribution form an IBVP. 

There are different approaches for modeling non-Darcy's phenomena [TTlfT^ll^ 
I18[ I20j . It can be derived from the more general Brinkman-Forchheimer's equation 
[18[ [7], or from mixture theory assuming certain relations between velocity field 
and "drag- like" forces due to fiuid to solid friction in the porous media |19j. It can 
be also derived using homogenization arguments 21 , or assuming some functional 
relation and then match the experimental data. In the current report, we just 
postulate a general constitutive equation relating the velocity vector field and the 
pressure gradient. We will introduce constraints on the momentum equation and 
on the fluid density. This will allow the reduction of the original system to a scalar 
quasi-linear parabolic equation for the pressure only. 

Hereafter, we use the following notations and basic definitions: 

• v{x,t) represents the velocity field; x is the spatial variable in M.'^,d = 2 
or 3; t denotes time; p{x,t) is the pressure distribution; y £ M.'^ are the 
variable vectors related to Vp; s, ^ represent scalar variables; 

• The notations C, Cq , Ci , C2 , . . . denote generic positive constants not de- 
pending on the solutions. 

Current studies of fiows in porous media widely use the three specific Forchheimer 
laws which we already mentioned (two-term, power, and three-term, respectively). 
Darcy and Forchheimer laws can be written in vector forms as follows: 

• Darcy's law: 

(1) av = — Vp, 

where a = — with k, in general, represents the permeability non- homogeneous 
k 

function depending on x subjected to the condition: fc^^ > k > ^2, 
1 > ^2 > 0. Here, the constant /i is the viscosity of the fluid. 

• The Forchheimer two-term law: 

(2) av + I3\v\v = -Vp, 

where /3 = , , F is Forchheimer's coefficient, $ is the porosity, and p is 

the density of the fluid. 

• The Forchheimer power law: 

(3) at. + c"|i;|"-iw = -Vp, 

where n is a real number belonging to the interval [1,2]. The strictly 
positive and bounded functions c and a are found empirically, or can be 
taken as c = (n — and a = a. By this way, n — 1 and n = 2 reduce 

the power law ([3]) to Darcy's law and to the Forchheimer two-term law, 
respectively. 

• The Forchheimer three-term law: 

(4) Av + B\v\v + C\vfv^-'^p. 
Here A, B, and C are empirical constants. 

We now introduce a general form for the Forchheimer equations. 

Definition II. 1 ((;- Forchheimer Equations). 

(5) g{x,\v\)v = ~\Jp, 
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here g{x, s) > for all s > 0. We will refer to ^ as g-Forchheimer (momentum) 
equations. 

Under isothermal condition the state equation relates the density p to the pres- 
sure p only, i.e. p = p{p). Therefore the equation of state has the form: 

(6) -:r"~' 

p dp K 

where 1/k is the compressibility of the fluid. For slightly compressible fluid, such 
as compressible liquid, the compressibility is independent of pressure and is very 
small, hence we obtain 



(7) P = Po exp 



P-Po 



K 



where po is the density at the reference pressure po (see [5] Sec. 2.3, and also [H] 
Sec. 3.4). Substituting Eq. ^ in the continuity equation 

(8) |-^^<'»). 

yields 

dp dp dp 

— — -pV • V - —V ■ Vp, 

dp dt dp 

(9) ^ = -kV - v-v - Vp. 

Since for most slightly compressible fluids in porous media the value of the constant 
K is large, following engineering tradition we drop the last term in ^ and study 
the reduced equation: 

(10) ^ = _kV-i;. 
^ ' dt 

II. 2. Boundary conditions. Let U C be a domain modeling the drainage 
area in the porous media (reservoir), bounded by two boundaries: the exterior 
boundary Fg, and the accessible boundary F^. 

The exterior boundary Fg models the geometrical limit of the well impact on 
the flow filtration and is often considered impermeable. This yields the boundary 
condition: 

(11) vN\r,=0, 

where N is the outward normal vector on the boundary F = F^ U Fg. Other types 
of boundary conditions on the exterior boundary are discussed in [4] . 

The accessible boundary Fi models the well and defines the regime of filtration 
inside the domain. On F^, consider a given rate of production Q{t), or a given 
pressure value p — ip{x,t), or a combination of both. It is very important from 
a practical point of view to build some "baseline" solutions capturing significant 
features of the well capacity and analyze the impact of the boundary conditions on 
these solutions. This analysis will be used to forecast the well performances and 
tune the model to the actual data. 

On the boundary F^, a "split" condition of the following type is of particular 
interest: 

(12) p^i;{x,t)=-f{t)+^{x), 
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where the time and space dependence of p are separated. This type of condition 
models wells which have conductivity much higher than the conductivity inside the 
reservoir. The limiting homogeneous case ipix) = const corresponds to the case of 
infinite conductivity on the well. 

In case the flow is controlled by a given production rate Q(t), the solution is not 
unique. 

Two important cases are: 

(a) pressure distribution of the form —At + f{x), and 

(b) constant total flux Q = const. 

The particular solutions of IB VP with boundary conditions (a) and (b) are "time- 
invariant" (see Section IIV|) and are used actively by engineers in their practical 
work. 

III. NON-LlNEAR DARCY EQUATION AND MONOTONICITY PROPERTIES 

In order to simplify the notation, we will further omit the a;-dependence of g in 
dS]). Thus, one has 

(13) g{\v\)^gix,\v\). 
From ([5]) one has 

(14) gi\v\)\v\ ^ ~\Vpl for s > 0. 

To make sure one can solve (jl4l) for \v\, we impose the following conditions: 
function g belongs to C([0, oo)) and C^((0,oo)), and satisfies 

(15) 5(0) > 0, and g'{s) > for all s > 0. 

Under this condition one has (55(5))' — sg'{s)+g{s) > g{0) > 0, for any positive 
value s. Therefore function sg(s) is monotone and one can find |?;| as a function of 
|Vp| 

(16) |^;| = G(|Vp|) 

Substituting equation ([T6l) into (O one obtains the following alternative form of 
the g-Forchheimer momentum equation ([5]): 

Definition III.l. (Non-linear Darcy Equation) [1] 

where the function K : [0,oo) — > [0,oo) is defined by 

(18) = 

Substituting ([TTl) for u into (fTO|) one derives the degenerate parabolic equation 
for the pressure: 

(19) ^^v(i^(|Vp|)Vp). 

where the compressibility constant k has been included in the non linear coefficient 
Ki\Vp\). 

Some properties of the function K can be found in [1]. It turns our that the 
function y — )■ K{\y\)y associated with the non-linear potential field on the RHS of 
equation (|17|) is monotonic. This monotonicity and related properties are crucial in 
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the study of the uniqueness and quahtative behavior of the the solutions of initial 
value problems (see e.g. [12]). 

We illustrate the function K for the particular case of two-term Forchheimer's 
equation. This is one of the few cases when the function K can be found explicitly. 

Example III. 2. For the Forchheimer two-term law (0), let g(s) = a + jSs, then 

one has sg{s) — /3s^ + as and s = G(^) = "^"^a" Thus 



213 

1 2 



We now introduce the notion of generalized polynomial with positive coefficients 
and positive exponents, abbreviated as GPPC, as it is a useful tool in this study. 

Definition III. 3. We say that a function g{s) is a GPPC if 

k 

(20) g{s) = aos"" + ais"i + ajs"^ + . . . + Uks"" = ^ ajs"^ , 

j=o 

where k > represents a natural number, = ao < ai < Q!2 < ■ • ■ < ftfc represent 
real values, and the coefficients ao,ai, . . . ,ak are real and positive. The largest 
exponent is the degree of g and is denoted by deg(g). 
Class (CPPC) is defined as the collection of all CPPC 

If the function g in Definition III. II belongs to class (GPPC) then we call it the 
g-Forchheimer polynomial. 

Lemma III. 4. Let g{s) be a function of class (CPPC) as in J^) . Then K{^) = 
Kg{^) is well-defined, is decreasing and satisfies 

(21) -^^<m)<^^.^^>o. 

where a = akjiak + 1) G [0, 1), and C'o and Ci are positive numbers depending on 
Gj 's and a J 's. 

The proof of the previous Lemma can be found in [1 and makes use of a condition 
that is automatically satisfied by a GPPC. 

Note that a = corresponds to the linear Darcy's case, while in the limiting case 
a — > 1 the largest exponent ak diverges. 

IV. PsEUDO Steady State Solutions and Productivity Index 

In engineering and physics, it is often essential to identify special time-dependent 
pressure distributions that generate flows which are time-invariant. In this section, 
we introduce the class of so-called pseudo-steady state (PSS) solutions which is 
used frequently by reservoir and hydraulic engineers to evaluate the "capacity" of 
the well (see. [HOHn] and references therein). 

Definition IV. 1. A solution p{x,t) of the equation I119\) in domain U, satisfying 
the Neumann condition on Pg is called a pseudo steady state (PSS) with respect to 
a constant A if 

(22) ^^W^ " " ~^ ^""^ 
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Equation ([T9| then implies 

(23) ^^^_A = V.(Al|Vp|)Vp). 

Using Green's formula and the Neumann boundary condition on Fg one derives 
A\U\ = -j {K{\\/p\)\/p) ■ Nda ^ j vNda = Q{t). 

Therefore, the total flux of a PSS solution is time-independent 

(24) Q{t) = A\U\=Q = const., for aU t. 

The PSS solutions inherit two important features, which we will explore further. 
On one hand, the total flux is defined by stationary equation (j23p and is given. On 
the other hand, the trace of the solution on the boundary is split a priori. Namely 
re-writing the PSS solution as 

(25) p{x,t) ^ -At + B + u{x), 

one has Vp = Vm, hence u and p satisfy the same boundary condition on Fg. On 
Fi, in general, we consider 

(26) p{x,t) = -At + B + Lp{x), on F,, 
where ^p{x) is given and the constant B is selected such that 

(27) j (pda = 0. 
Therefore u{x) satisfies 

(28) - A = V ■ K{\Vu\)Vu, 
, , du 

(30) u — ip on Fj. 

Of particular interest is the case 1^9(2;) = on F^. From a physical point of view, 
it relates to the constraint that conductivity inside well is non-comparably higher 
than in the porous media. 

We call u{x) the profile of PSS corresponding to A and the boundary profile 

ip{x). 



Remark IV. 2. Note that for a PSS as in i25\) . we have the quantity 

Q{t) _ Q 



J{t) 



|FT lu P(^^ ^)^^ " ir-T /r, Pi^^ W\ lu u{x)dx ■ 



It represents the production rate versus the pressure drawdown (the difference be- 
tween averages in the domain and on the boundary Ti), and is independent of time. 
This quantity is called Productivity Index, and it is widely used by engineers to test 
the performances of a well/reservoir system (see [4l l3l [15] j. 
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V. Geometric Interpretation 

In the PDE literature equation (see Serrin G., 1967, Gilbarg, D. & Trudinger, 
N.1977, Evans L. 1999) 

(31) div (jl + \Wuf^ 

is referred as homogeneous CMC equation were u is a graph defined in the domain 
U (Z ■ We win also refer to it as a CMC graph equation. We noted that this 
equation looks somewhat similar to 

(32) div{K{\\Iu\)Vu) = -A, 

which has been previously introduced in defining the basic PSS profile for the non- 
linear Forchheimer equation. 

It is not difficult to see that the non linear term in (|3T|) is "about" (1 + |Vm|)~i, 
which means that there exist constants Cq and Ci such that 

(33) T^.<-, '—T7,< 



1 + |Vm| - / , ,2\i/2 - l + |Vu| 



(i + |v^l') 



We also showed in Lemma [111.41 that the non linear term X(|Vp|) in (15^ is, in the 
GPPC case, "about" (1 + |Vu|°)-i. Here a = ak/{ak + 1) < 1, where > 
is the degree, deg{g), of the specific GPPC polynomial. Then in the limiting case 
a — 1, we can expect the two non- linear coefficients in Eqs. (|311) and (|32p to have 
the same structure. Therefore one cannot expect the PSS Forchheimer equation to 
"survive" in limiting case a = 1. At the same time it is worth mentioning that, for 
the case a < 1, the weak solution of the PPS Forchheimer equation is unique (see 
[T]) and exists in the corresponding Sobolev space W^'^~'^. This result was treated 
in detail in [T]. 

In the next section we will introduce a few basic geometric notions and defi- 
nitions, in order to show the robust link between these two objects. The actual 
relationship between the two equations is far from being straight-forward, in spite 
of their similarity. In particular, we will show that the pressure function can be 
interpreted as a CMC graph under some constraints and appropriate geometric 
transformations . 



VI. The mean curvature equation for a graph. 

Any map r : D C whose differential map dr is injective is called 

an immersion (surface immersion) in R^. Equivalently, the map r represents an 
immersion if its Jacobian has rank 2, or all points are regular. If an immersion r 
is 1-1, it is sometimes called an embedding. Any immersed surface can be endowed 
with a general Riemannian metric g — g{x,y) (|22j. page 418). In the particular 
case when the metric g is defined using the usual velocity vector fields, we will call 
it naturally induced metric (i.e., naturally induced by the immersion). In the most 
usual notation, (M,g) represents a Riemannian manifold of Riemannian metric g, 
and M = r{D). 

Consider a smooth surface that can represented as a graph z — u{x,y) of an open 
domain D <Z R? . This surface is parameterized via the map r : D <Z R? ^ R^ ^ 

(34) r{x,y) ^ {x,y,u{x,y)) . 
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Definition VI. 1. (Natural metric) We will call naturally induced metric the fol- 
lowing quadratic differential form: 

(35) dr^[x,y) ^\r^\^dx^ + 2 <r^,ry> dx dy + \ry\^dy^ , 
which can he rewritten as 

(36) dr'^{x, y) = {1 + u^) dx^ + 2uxUy dx dy + {\ + Uy) dy^ , 
where the coefficients 

511 = 1 + , 512 = Uy 522 = 1+ 

represent the entries of the corresponding matrix g. 

Definition VI. 2. (Gauss Map) We will call Gauss map the usual unit normal 
vector field N : D ^ defined as: 

rr, X r„ —Mr i — u,, i + k 

II ^ II ^ul + ul + l 

Definition VI. 3. (Second Fundamental form) The second fundamental form is 
defined by 

(38) da'^{x, y) — hu dx"^ + 2 hi2 dx dy + /i22 dy'^ , 
with the following coefficients 

(39) /ill :=< 7V,r^:r >, hi2 ■.^< N,r^y > , ■■^< N,ryy >, 
and so it can be expressed as 

(40) da^ix, y) = , dx^ + 2 dx dy + dy^ . 

Note that once the local coordinates are given, one can identify the first and 
second fundamental forms with their corresponding 2x2 matrices, g and ft., re- 
spectively. The matrix operator S = g~^h can be viewed as a linear operator from 
the tangent plane to the surface at a given point, to the same tangent plane. S is 
usually called the Weingarten map or shape operator (see [22], vol. II). Note that 
5'-' • hij represents its trace. The following result is a classical result of differential 
geometry, which can be found in any text-book (e.g., [22], vol. II), and whose proof 
is elementary. 

Proposition VI. 4. For any first and second fundamental forms defined for an 
immersion r ~ r{x,y), the following formula is satisfied: 

(41) 5^^ • h,, = 2H, 

where g and h are matrices corresponding to the first and second fundamental form, 
and H represents the mean curvature, defined as the arithmetic mean of the prin- 
cipal curvatures of the immersion r. 



This equation is frequently referred to as mean curvature equation. 
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Definition VI. 5. (Laplace- Beltrami operator) The Laplace- Beltrami operator AgU 
corresponding to the graph z — u{x, y) and the metric defined in eq \3(j]) is defined 
as (see 



(42) Ag u = ^_ • g 



^ \/dcig ^ dx^ dx^ 
It is worthwhile noting that the Laplace-Beltrami operator is frequently defined 

without the factor ; on the other hand, this factor plays an important role 

Vdetg 

in our work. Historically speaking, its original definition is the same as ours. Also, 
it is important to clarify that the Laplace Beltrami operator of a general surface 
immersion r is defined component-wise and represents a vector (r) , which in our 
case it reduces to its last component, lS.g(u). 

[9] was the first well-known reference to make the observation that for surfaces 
immersed in the operator g"^^ ■ hij acting at each point coincides with the Laplace- 
Beltrami operator. In view of this observation and the previous definition, the mean 
curvature equation (13) can be rewritten as 



Ag u = 2H. 

A immersion r : D C ^ with vanishing Laplace-Beltrami operator is 
said to be harmonic in a generalized (Riemannian) sense. 

Remark that for the case when the surface metric is the flat Euclidean one 
g represents the identity matrix and the Laplace-Beltrami operator becomes the 
usual Laplace operator A — dxx + dyy . 

A well-known result of geometric surface theory stated that an immersion r as 
above is harmonic if and only if it is minimal [5]. This result agrees with our 
previous Proposition and Definition. 

The mean curvature equation will be referred to as constant mean curvature 
equation for the case when H is constant. The case oi H = const ^ represents 
the case of CMC surfaces (to be distinguished from the minimal case H = 0). 

We showed in P] that one can independently modify the coordinate functions of 
the velocity vector, in a way that links the CMC equation to the Forchheinier-type 
equation. 

We started from the partial velocities rx and ry at the point P{x,y) of the initial 
graph r{x^y) = {x,y,u{x,y)), and we applied the following transformation: 

(43) rx ^ f-x -.^ {x,0,lJ-{x,y)ux(x,y)) 

(44) ry ry (0, x, ^x, y)uy{x, y)) , 

where x is a scaling constant and ^{x, y) is a smooth function. 

Note that the notations fx , fy do not a priori assume the existence of an immer- 
sion f whose partial velocity fields are written in this form. We have studied the 
conditions in which there exists such a parametric surface. This represents an easy 
case of the Frobenius theorem, and the existence condition of such an immersion 
f reduces to the compatibility condition (fx)y = {fy)x being satisfied for the vec- 
tors defined above. For details, see Example 1.2.3, from [16], based on successively 
applying Picard's Theorem in the x and y-directions, respectively. Assuming that 
this compatibility condition is satisfied. Example 1.2.3, from [16], shows that for 
any fixed initial position Pq at the origin (0,0) (or another fixed base point), there 
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exists a unique solution r = {x , y , u{x , y)) , as an immersion whose partial velocity 
vectors are fx and fy. 

The change u{x, y) — u{x^ y) represents a smooth deformation of the graph along 
the z-axis, while the change {x,y) — >■ with x = ^-nd y = XV represents 

the rescaling of the graph domain from D to D. This non-trivial transformation 
modifies the Gauss map (the tangent plane), as well as the shape operator. 

Proposition VI. 6. Consider a smooth graph in that is viewed as an immersion 
from an open simply connected planar domain D into the Euclidean space via 
r{x,y) = {x,y,u{x,y)). 

At every point on M — r{D) consider the modified velocity vectors defined as: 

fx = (x, 0, ^l{x, y)ux) and fy = (0, x, /i(a;, y)uy), 

where x constant and fi{x, y) is an arbitrary smooth function such that fx 
and fy are non-vanishing and linearly independent and such that the compatibility 
condition (fx)y — {fy)x is satisfied (or equivalently, fixUy — fiyUx). 

Consider a fixed initial position Pq at the origin (0,0) (or base point). Let 
f{x,y) = {x,y,u{x,y)) represent the uniquely determined integral surface having fx 
and fy as partial velocity vectors, such that f(0,0) = Pq. These vectors naturally 
induce the first fundamental form g . The corresponding acceleration vectors together 
with the Gauss map N{x,y) naturally induce the second fundamental form h. 

Then the corresponding Laplace- Beltrami operator can be interpreted in terms of 
the trace of the shape operator, that is, 



Remark VI. 7. Note that f actually represents a family of immersions of parameter 
/i, a real valued smooth function of two variables. The initial immersion r belongs 
to this family, corresponding to the case X = M — 1- call z = u(x,y) 

generalized graph, and g its naturally induced metric. 

Example 1.2.3, from [16j, based on successively applying Picard's Theorem in 
the X and y-directions, provides an explicit solution for f above. Assuming that the 
compatibility condition is satisfied, we denote by f = {x,y,u(x,y)) an immersion 
whose partial velocity vectors are fx and fy. Some computational details can be 
found in [2]. We therein collect the main information on the first and second 
fundamental forms of the generalized graph, namely: 

We first derive the expressions of the data corresponding to the immersion f. 
The induced metric is given by 

(46) df'^{x, y) = (x^ + t'-'^u^x) + 2y.^Ux Uy dxdy + (x^ + f-i^Uy) dy^ . 
The unitary normal vector field is given by 



(45) 



(47) 



N = 




liUxi- lJ.Uyj-\-x'^ 
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The coefficients of tlie second fundamental form, as entries of the corresponding 
matrix h, are: 

(48) hn^ Xif^u^. + u^t^..) 

(49) h,2 = 

(50) h22 = 



X(M"w + %Mj;) 



The corresponding Laplace-Beltrami operator is given by 
(51) 

_ (X^ + A^^ "y) Mg^o: - 2// Ua; My + (x^ + "x) "to , xi^x fJ-x + Uy ^ly) 

L\a - — 



X (X^ + («? + tig))'/' ' (x2 + (7X2 + «2))3/2 

which can be rewritten as 



(52) A,.V. "J" Uv^U^|.2g. 



Here V- and V are the divergence and the gradient operator in the new reference 
system (i, y). 

Remark VI. 8. In case H is constant, u is a CMC graph with respect to the domain 
D, provided a solution to L5S^) exists. 

The following theorem is an immediate consequence (see [2j for details of the 
proof): 

Theorem VI. 9. Consider a smooth graph u{x,y) in the domain D[x,y), solution 
of 

(53) V ■ {K{\Vu\)Vu) = -A, 
and such that 

(54) {Ux Uxy ~t~ Uy Uyy'^ Ux — {^Ux Uxx ^yx^ ^y ■ 

Let D{x,y) be the scaled domain obtained by the conformal mapping x — xx and 
y ^ XD- ^(^j y) ^6 stretched graph, which is parameterized as the immersion 
f with partial velocities (fx,fy) given by J^'^K^-^l i and 

(55) M(x,|V.|) = ^=i£M_, 
^ ' v/l-x'^(|Vw|)2|Vu|2 
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(56) X < Xma: 



max ^ max 

Then u{x^y) is solution of the corresponding CMC equation 
(57) V- I I ^A. 



1 + IVu 



r,|2 
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Remark VI. 10. Condition J56)) assures fi to be a real smooth valued function. 
Condition (j54p assures the compatibility condition 

which are necessary and sufficient conditions for u to exist. Reformulated, (j54p 
states that each level curve of the graph of u (at z = ci) represents a level curve of 
the graph of |Vu| (at z — C2)- 

Replacing the generic function fi with ([55|) specified in the statement of this 
theorem into equation (1521) immediately gives 



(58) V ■ \ ^^^=\ ^-y ■ {Ki\Vu\)Vu) 

This equality practically maps solutions of the Forchheimer equation (I53p to 
solutions of the CMC equation ([5^ through an explicit transformation, under a 
certain natural assumption. 

Proposition VI. 11. Consider a smooth graph u{x,y) solution of the PSS Forch- 
heimer equation I15S\) and satisfying ^54^ . Consider the associated CMC graph 
u{x, y) solution of |57p with /i given by ( f?5|) and x < Xmax- Let rj = |Vu|, £, — iV-uj, 
and T = . ^ , then 

(59) n = 9{\v{T,x)\)HT,x)\, 
where 

(60) |«(r,x)| = -. 

Here v is the velocity of the fluid flow in the porous media. As a consequence 

(61) = -V-j- 

(62) uy^^v^ 

Proof. The proof is based on the vector identity 

^(|Vu|)Vu = Vu, 
which implies the scalar equality 

(M|Vu|)|V«|)2 = |VSp. 
Taking into account equation (j55p for fi, and substituting rj for |Vm| and ^ for 



|Vu|, we obtain 

(63) TTTTT-T-^^^ • 



Let us consider r = xK{r])r], which is positive. 
Substituting and solving for r in we obtain 

(64) t(0 = , ^ . 
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Recollecting Eq. (fT7|) for the velocity field, it follows 

(65) \v{T,x)\^Kir^)r^^-. 

X 

Now by using the g-Forchheimer Eq. (III.ip . it follows Eq ((59)) . Clearly, the corre- 
spondence between ^ and rj is one-to-one. □ 

VI. 1. Application to reservoir engineering. The following proposition was 
proved in j^. We recall it here only because the proposed solution strategy for 
the evaluation of the productivity index is among the direct applications of the 
results hereby presented. 

Proposition VI. 12. For the GPPC case the PSS Productivity Index can he com- 
puted as 



(66) P/(a„a„|«|) 



Jgi\v\)\v^dx JELoaM^'^^d^' 

u u 



Here PI is time invariant, and it depends explicitely on ai,ai , \v\ and domain U. 

This result combined with previous theorem (in particularly with Eq. ()60p for 
the velocity) expresses the fact that the Productivity Index of the well can be eval- 
uated in the following way, provided that condition (j54p is verified. 

Solution strategy for the evaluation of the Productivity Index: 

(1) The factor x is selected to generate the scaled domain D. 

(2) The CMC equation ([57]) is solved for u{x, y) on the domain D for Q, with 
the appropriate boundary conditions. 

(3) The gradient norm ^ = iVwl is evaluated. 

(4) The coefhcients a^, a^, i = 0, - ■ ■ ,k of the GPPC polynomial are selected. 

(5) The norm of the velocity \v\ is explicitely evaluated trough Eq. (|60|) . 

(6) The Productivity Index PI{ai, ai, \v\) is evaluated by using formula (|66p . 



Remark VI. 13. Points 1 and 2 do not depend by the choice of ai and ai which 
means that the evaluation of in 3 is and ai independent. We need to solve just 
one EVP. The Productivity Index can be evaluated a posteriori for any choice Oi 
and ai. 

VII. Acknowledgments 

The authors would like to thank Luan Hoang for his stimulating discussions. 
This research is supported by NSF Grant No. DMS-0908177. 

References 

[1] E. Aulisa, L. Bloshanskaya, L. Hoang, and A.Ibragimov, Analysis of generalized Forchheimer 
flows of compressible fluids in porous media. J. Math. Phys. Vol. 50, Available online (2009). 

[2] E. Aulisa, A. Ibragimov, M. Toda, Geometric Frame- Work for Modeling Non-Linear Flows 
in Porous Media and Its Applications in Engineering, Journal of Non-linear Analysis - 
Real World Applications, vol. 10 (2009); currently in press, avail, at Science Direct, 
doi:10.1016/j.nonrwa.2009.03.028 

[3] E. Aulisa, A. L Ibragimov, P. P. Valko, J. R. Walton, Mathematical Frame-Work For Produc- 
tivity Index of The Well for Fast Forchheimer (non-Darcy) Flow in Porous Media. Mathe- 
matical Models and Methods in Applied Sciences, DOL10.1142/S0218202509003772, in press. 



GEOMETRIC METHODS IN THE ANALYSIS OF NON-LINEAR FLOWS IN POROUS MEDIA5 

[4] E. Aulisa, A. I. Ibragimov, P. P. Valko, J. R. Walton, A new method for evaluation the 
productivity index for non-linear flow, SPE Journal, DOI 10.2118/108984-PA, in press (2009). 

[5] J. Bear, Dynamics of Fluids in Porous Media. Dover Publications, Inc., New York, 1972. 

[6] L. Bloshanskaya, E. Aulisa and A. Ibragimov, Mathematical model of well productivity index 
for generalized Forchheimer flows and application. Proceedings of SIAM Conference Mathe- 
matics for Industry: Challenges and Frontiers (MI09). 

[7] J. Chadam, Y. Qin, Spatial decay estimates for flow in a porous medium, SIAM J. Math. 
Anal., Vol. 28, No. 4, 808-830 (1997). 

[8] L. P. Dake, Fundamental in reservoir engineering, Elsevier, Amsterdam, (1978). 

[9] J. Eells, L. Lemaire, Two reports on harmonic maps, World Scientific Pub Co Inc, July 1995, 
ISBN-10: 9810214669. First report. 
[10] J. Eells, A. Ratto, Harmonic maps and minimal immersions with symmetries. Methods of 
Ordinary Differential Equations Applied to Elliptic Variational Problems, Annals of Mathe- 
matics - Studies - Princeton University Press, Study 130, 1993. 
[11] E. Ewing, R. Lazarov, S. Lyons, D. Papavassiliou, Numerical well model for non Darcy flow, 

Comp. Geosciences, 3, 3-4, 185-204 (1999). 
[12] L. C. Evans, Partial Differential Equations. American Mathematical Society, Providence, 
(1998). 

[13] P. Forchheimer, Wasserbewegung durch Boden Zeit, Ver. Deut. Ing. 45, (1901). 

[14] M. Muskat, The flow of homogeneous fluids through porous media, McGraw-Hill Book Com- 
pany, Inc., New York and London, (1937). 

[15] A. I. Ibragimov, D. Klialmanova, P. P. Valko, J. R. Walton, On a mathematical model of the 
productivity index of a well from reservoir engineering, SIAM J. Appl. Math., 65, 1952—1980 
(2005). 

[16] T. A. Ivey, J. M. Landsberg, Cartan for Beginners, Differential Geometry of Moving Frames 
and Exterior Differential Systems, Graduate Studies in Mathematics, AMS, vol. 61, 2003. 

[17] J.L. Miskimins, H.D. Lopez-Hernandez, R.D. Barree, Non-Darcy Flow in Hydraulic Frac- 
tures, JPT, March, pp. 57-59, 2006 

[18] L. E. Payne, B. Strauglian, Convergence and Continuous Dependence for the Brinkman- 
Forchheimer Equations, Studies in Applied Mathematics, 102, 419-439 (1999). 

[19] K. K. Rajagopal, On a hierarchy of approximate models for flows of incompressible fluids 
through porous solids, Mathematical Models and Methods in Applied Sciences, Vol. 17, No. 
2, 215-252 (2007). 

[20] K. R. Rajagopal, L. Tao, Mechanics of mixtures, World scientific, (1995). 

[21] E. Sanchez-Palencia, Non-Homogeneous Media and Vibration Theory, Lecture Notes in 

Physics, Springer- Verlag, (1980). 
[22] M. Spivak, Differential geometry, vol. I and II, 2nd Edition, Publish or Perish, 1977. 
[23] C.A.P. Tavera, H. Kazcmi, E. Ozkan, Combine effect of Non-Darcy Flow and Formation 

Damage on Gas Well Performance of Dual-Porosity and Dual Permeability Reservoirs, SPE- 

90623, 2004. 

[24] S. Whitaker, The Forchheimer Equation: A Theoretical Development, Transport in Porous 
Media 25: 27-61 (1996). 

Department of Mathematics and Statistics, Texas Tech University, Box 41042, Lub- 
bock, TX 79409-1042, U. S. A. 

E-mail address: eugenio.aulisaQttu.edu 
E-mail address: akif.ibraguimov@ttu.edu 
E-mail address: magda.todaQttu.edu 



